********************************************************************************
***************** Alternative behavioral factors affecting the wedge

clear all
do "$maindir/Code/bootstrapres.do"

* CBA information based on marginal social costs
use "$maindir/Results/Model_Outputs/ihwap_wedge.dta"

*** prism-related variables
gen rsq_pre = max if treated==0
bysort Household : egen rsquare_pre = max(rsq_pre)
gen rsq_post = max if treated==1
bysort Household : egen rsquare_post = max(rsq_post)

gen hd_pre = HDDbest if treated==0
bysort Household : egen hdd_pre = max(hd_pre)
gen hd_post = HDDbest if treated==1
bysort Household : egen hdd_post = max(hd_post)

* prism sample restriction
gen prismrest = 0 if rsquare_pre!=. & rsquare_post!=. & hdd_pre!=. & hdd_post!=.
replace prismrest = 1 if hdd_pre>=43 & hdd_pre<73 & hdd_post>=43 & hdd_post<73 & ///
		rsquare_pre>=0.85 & rsquare_pre!=. & rsquare_post>=0.85 & rsquare_post!=.

** merge with info about furnace replacements
merge m:1 Household using "$maindir/Data/furn_replace.dta"
drop if _merge==2
drop _merge


***** some extra restrictions on number of homes served by each contractor
gen aux=1
sort Household meterend
bysort Household : gen homeobs = sum(aux) if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs = . if homeobs!=1

sort Household meterend
bysort Household : gen homeobs_post = sum(aux) if treated==1 & twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009
replace homeobs_post = . if homeobs_post!=1

** fixing builddate to avoid omitted variables
replace builddate = 11 if tab_builddate11==1 | tab_builddate12==1 | tab_builddate13==1

** fixing sqfeet to avoid omitted variables
drop roundsqfeet
gen roundsqfeet = .
	forval i = 0(500)2500 {
		replace roundsqfeet = `i' if sqfeet>`i' & sqfeet<=`i'+500
	}
replace roundsqfeet = 3000 if sqfeet>3000 & sqfeet!=.

** fixing some spending variables to avoid omitted bins
drop binAirCon
gen binAirCon = .
replace binAirCon=0 if Real_tot_actAirCon==0
replace binAirCon=1 if Real_tot_actAirCon>0 & Real_tot_actAirCon!=.

replace binAttic = 2700 if Real_tot_actAttic>2400 & Real_tot_actAttic!=.

** demographic variables
*tabulate group_demographics, gen(tab_demog)
*drop tab_demog14
gen simpocc = noccupants
replace simpocc = 5 if noccupants>=5 & noccupants!=.

gen binAge = .
forval i = 20(10)80 {
replace binAge = `i' if Age>`i' & Age<=`i'+10
}

gen binIncome = .
forval i = 5000(5000)35000 {
replace binIncome = `i' if Real_income>`i' & Real_income<=`i'+5000
}
replace binIncome = 1 if Real_income<=5000
replace binIncome = 40000 if Real_income>40000 & Real_income!=.


*** housing structure
gen binSqft = .
forval i = 600(300)2700 {
replace binSqft = `i' if sqfeet>`i' & sqfeet<=`i'+300
}
replace binSqft = 1 if sqfeet<=600
replace binSqft = 3000 if sqfeet>3000 & sqfeet!=.

drop BlowerPreBins
gen BlowerPreBins = .
forval i = 2000(1000)9000 {
replace BlowerPreBins = `i' if Blower_Pre>`i' & Blower_Pre<=`i'+1000
}
replace BlowerPreBins = 1 if Blower_Pre<=2000
replace BlowerPreBins = 10000 if Blower_Pre>10000 & Blower_Pre!=.

*** energy prices
gen binelec = .
local counter = 2
forval i = 10.5(0.5)12.5 {
replace binelec = `counter' if elec_prices>`i' & elec_prices<=`i'+0.5
local counter = `counter' + 1
}
replace binelec = 1 if elec_prices<10.5
replace binelec = `counter' if elec_prices>13 & elec_prices!=.

gen bingas = .
local counter = 2
forval i = 7(1)16 {
replace bingas = `counter' if gas_prices>`i' & gas_prices<=`i'+0.5
local counter = `counter' + 1
}
replace bingas = 1 if gas_prices<7
replace bingas = `counter' if gas_prices>17 & gas_prices!=.

* transforming decimals in percent
gen savings_gap_pct = savings_gap*100

forval i = 1/200 {
gen savings_gap_pct`i' = 100*savings_gap`i'
}

** counting number of homes served by contractor
merge m:1 Household using "$maindir/Data/ihwap_state.dta", keepusing(arch_contID) nogen keep(3)
sort arch_contID
by arch_contID : egen totjobs = total(aux) if ww_treat==1
by arch_contID : egen totaljobs = max(totjobs)
drop totjobs
replace totaljobs=. if arch_contID==.

*br if totaljobs>700 & totaljobs!=.
* Contractor with ID 236 to be ommitted, being one of the biggest contractors in the sample

** main heat pre-treatment
gen MainHeatkBTU = MainHeatBTU/1000

gen binheatkbtu = .
forval i = 60(10)120 {
replace binheatkbtu = `i' if MainHeatkBTU>`i' & MainHeatkBTU<=`i'+10
}
replace binheatkbtu = 1 if MainHeatkBTU<=60
replace binheatkbtu = 130 if MainHeatkBTU>130 & MainHeatkBTU!=.

gen heatworks = 0 if MainHeatOperational=="False"
replace heatworks = 1 if MainHeatOperational=="True"


* fixing variables that have a zero indicator (for later interacting)
foreach x in nstories binAirCon white black hispanic otherrace haselderly hasminor heatworks {
replace `x' = `x'+1
}
replace roundAtticR = 1 if roundAtticR==0

gen heatnotwork = 1 if MainHeatOperational=="False"
replace heatnotwork = 0 if MainHeatOperational=="True"

replace binFurnace = 1 if binFurnace==0

gen failprism = 0 if prismrest==1
replace failprism = 1 if prismrest==0

* variables for projected savings
gen ExistingConsumption_mmbtu = ExistingConsumption/1000000
gen ProjectedConsumption_mmbtu = ProjectedConsumption/1000000


***** climate zones
merge m:1 Household using "$maindir/Data/wx_weather.dta", keepusing(region)
drop if _merge==2
drop _merge
egen regionID = group(region)


********************************************************************************
*** correlate PRISM HDD and r-squared
twoway (scatter hdd_pre rsquare_pre , msize(vsmall) mcolor(black)) , ///
	ytitle("PRISM Optimal HDD (Pre-Treatment)") ///
	xtitle("PRISM R-Squared (Pre-Treatment)") ///
	graphregion(color(white)) bgcolor(white) ///
	xline(0.85, lpattern(dash))	yline(43, lpattern(dash)) yline(72, lpattern(dash)) ///
	legend(off)
graph export "$maindir/Results/Figures/corr_hdd_rsq_pre.png", replace width(2500)

twoway (scatter hdd_post rsquare_post , msize(vsmall) mcolor(black)) , ///
	ytitle("PRISM Optimal HDD (Post-Treatment)") ///
	xtitle("PRISM R-Squared (Post-Treatment)") ///
	graphregion(color(white)) bgcolor(white) ///
	xline(0.85, lpattern(dash))	yline(43, lpattern(dash)) yline(72, lpattern(dash)) ///
	legend(off)
graph export "$maindir/Results/Figures/corr_hdd_rsq_post.png", replace width(2500)



********************************************************************************
******* how wedge is affected by other behavioral factors

**** specification omitting ZERO spending
reg savings_gap_pct 1.heatnotwork 1.failprism ib1.binFurnace ib1.binAirSeal ib1.binAttic ib1.binWallIns ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		ib1.white ib1.black ib1.hispanic ib1.otherrace ib1.haselderly ib1.hasminor ///
		ib1500.binSqft ib1.nstories ib1.roundAtticR ib1.builddate ///
		ib2011.ProgramYear ib236.arch_contID ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009 & treated==1 & totaljobs>=20 & totaljobs!=. ///
		, vce(cluster Household) 	
est sto gapreg
matrix gapbetas = e(b)
matrix gapvars = e(V)
sca nobs = e(N)
sca rsq = e(r2)
sca rss_full = e(rss)


local dim `= rowsof(gapbetas)',`=colsof(gapbetas)'
matrix gap_coefs_boot2 = J(`dim',.)
local colnames : colnames gapbetas
mat colnames gap_coefs_boot2=`colnames'
mat rownames gap_coefs_boot2="y1"
forval i = 1/200 {
	di "starting bootstrap iteration `i'"
	qui : reg savings_gap_pct`i' 1.heatnotwork 1.failprism ib1.binFurnace ib1.binAirSeal ib1.binAttic ib1.binWallIns ///
		ib1.binBaseload ib1.binDoor ib1.binWtHtr ib1.binHealSfty ///
		ib1.binAirCon ib1.binFoundation ib1.binGeneral ib1.binWindow ///
		ib3000.BlowerPreBins ib70.binheatkbtu ///
		ib2.simpocc ib50.binAge ib15000.binIncome ///
		ib1.white ib1.black ib1.hispanic ib1.otherrace ib1.haselderly ib1.hasminor ///
		ib1500.binSqft ib1.nstories ib1.roundAtticR ib1.builddate ///
		ib2011.ProgramYear  ib236.arch_contID  ///
		c.tmin c.tmin#c.tmin c.tmax c.tmax#c.tmax c.precip c.precip#c.precip ///
		i.end_month elec_prices gas_prices ///
		[fweight = weight_boots`i'] ///
		if twoyears_prepost==1 & meterdays_monthly!=. & ProgramYear>=2009 & treated==1 & totaljobs>=20 & totaljobs!=.

	matrix tempmat = e(b)
	mat_rapp  gap_coefs_boot2 : gap_coefs_boot2 tempmat, miss(.)
}


****** output coefficients with bootstrapped stanard errors
matrix gapcoefs = gapbetas
matrix bootstrapcoefs2 = gap_coefs_boot2[2..., 1..332]
mata : st_matrix("bootstrapcoefs2", editvalue(st_matrix("bootstrapcoefs2"), 0, .))
mata : st_matrix("gap_var_boot2", mm_colvar(st_matrix("bootstrapcoefs2")))
matrix gap_var_boot2 = diag(gap_var_boot2)
mata : st_matrix("gap_var_boot2", editvalue(st_matrix("gap_var_boot2"), ., 0))
local colnames : colnames gapbetas
mat colnames gap_var_boot2=`colnames'
mat rownames gap_var_boot2=`colnames'
bootstrapres gapcoefs gap_var_boot2 nobs
est sto decompose_gap_behavior


esttab decompose_gap_behavior ///
				using "$maindir/Results/Tables/gapreg_behavior.tex", ///
				replace b(3) se(3) nostar nonumbers ///
				nonotes stats(N, fmt(%9.0gc) labels("Observations")) ///
				keep(1.heatnotwork 1.failprism)


